Crosstalk between 5-methylcytosine and N6-methyladenosine machinery defines disease progression, therapeutic response and pharmacogenomic landscape in hepatocellular carcinoma

Background Accumulated evidence highlights the significance of the crosstalk between epigenetic and epitranscriptomic mechanisms, notably 5-methylcytosine (5mC) and N6-methyladenosine (m6A). Herein, we conducted a widespread analysis regarding the crosstalk between 5mC and m6A regulators in hepatocellular carcinoma (HCC). Methods Pan-cancer genomic analysis of the crosstalk between 5mC and m6A regulators was presented at transcriptomic, genomic, epigenetic, and other multi-omics levels. Hub 5mC and m6A regulators were summarized to define an epigenetic and epitranscriptomic module eigengene (EME), which reflected both the pre- and post-transcriptional modifications. Results 5mC and m6A regulators interacted with one another at the multi-omic levels across pan-cancer, including HCC. The EME scoring system enabled to greatly optimize risk stratification and accurately predict HCC patients’ clinical outcomes and progression. Additionally, the EME accurately predicted the responses to mainstream therapies (TACE and sorafenib) and immunotherapy as well as hyper-progression. In vitro, 5mC and m6A regulators cooperatively weakened apoptosis and facilitated proliferation, DNA damage repair, G2/M arrest, migration, invasion and epithelial-to-mesenchymal transition (EMT) in HCC cells. The EME scoring system was remarkably linked to potential extrinsic and intrinsic immune escape mechanisms, and the high EME might contribute to a reduced copy number gain/loss frequency. Finally, we determined potential therapeutic compounds and druggable targets (TUBB1 and P2RY4) for HCC patients with high EME. Conclusions Our findings suggest that HCC may result from a unique synergistic combination of 5mC-epigenetic mechanism mixed with m6A-epitranscriptomic mechanism, and their crosstalk defines therapeutic response and pharmacogenomic landscape. Supplementary Information The online version contains supplementary material available at 10.1186/s12943-022-01706-6.

extrinsic and intrinsic immune escape mechanisms, and the high EME might contribute to a reduced copy number gain/loss frequency. Finally, we determined potential therapeutic compounds and druggable targets (TUBB1 and P2RY4) for HCC patients with high EME. Conclusions Our findings suggest that HCC may result from a unique synergistic combination of 5mC-epigenetic mechanism mixed with m 6 A-epitranscriptomic mechanism, and their crosstalk defines therapeutic response and pharmacogenomic landscape.
Keywords Hepatocellular carcinoma, 5-methylcytosine, N 6 -methyladenosine, Therapeutic response, Pharmacogenomic landscape, Multi-omics Background Liver cancer, of which 90% are hepatocellular carcinoma (HCC), is the seventh most frequent cancer globally, with 905,677 new diagnosed cases (4.7%) and 830,180 (8.3%) new death cases globally according to the GLOBOCAN 2020 statistics [1]. Most HCC cases arise in the context of chronic liver disease or cirrhosis, the most common cause of which is non-alcoholic fatty liver disease, alcohol-related liver disease, and hepatitis B virus (HBV)/hepatitis C virus (HCV) infection [2]. HCC management is defined by the Barcelona Clinic Liver Cancer (BCLC) staging system that bases on tumor burden, liver function, physical status, etc. [3]. Hepatectomy, ablation and liver transplantation are recommended for patients at the very early stage or early stage, with transarterial chemoembolization (TACE) for those at the intermediate stage, and systemic therapy for those at the advanced stage with adequate liver function and good performance status [3]. Despite the promising treatment options, prognostic outcomes of HCC are still bleak due to the high risk of recurrence and metastasis, with a 5-year survival rate of no more than 12% for advanced HCC [4]. Sorafenib, a tyrosine kinase inhibitor with angiogenic and proliferative effects, is widely applied in the systemic therapy of advanced HCC. However, many HCC patients are not responsive to sorafenib or become resistance within 6 months. Immuno-oncology has become a paradigm shift for the treatment of human cancers, including HCC. The IMbrave150 clinical trial showed that PD-L1 inhibitor (atezolizumab) in combination with VEGF inhibitor (bevacizumab) presented the superior progression-free survival (PFS) and overall survival (OS) in comparison to sorafenib among unresectable HCC patients [5]. Although the combined immunotherapy has been approved as the frontline therapeutic option for advanced HCC patients, the response rate remains ~ 30%. Due to the wide heterogeneity of risk factors and pathogenesis of HCC, established strategies for prediction and prognostication remain limited [6]. Hence, it is critical to ascertain new methods to improve the early diagnosis of HCC and to predict treatment response and survival in patients with established HCC.
Cellular DNA and RNA undergo various forms of methylation. Methylation of DNA and RNA, especially 5-methylcytosine (5mC) and N 6 -methyladenosine (m 6 A) modifications exert key roles in a variety of biological processes [7]. DNA methylation is a well-known and critical epigenetic modification. Studies have revealed that 5mC is the most common type of DNA modification in eukaryotes [8]. m 6 A is the most abundant form of mRNA modification in eukaryotes, which regulates several processes of mRNA metabolism, especially mRNA translation and degradation [9]. Accordingly, several 5mC and m 6 A machines have been well identified. Based on the diverse functions of 5mC and m 6 A machines, 5mC and m 6 A modifications have been shown to impact many fundamental biological processes, including HCC [10,11]. Rapidly accumulating evidence demonstrates the significant crosstalk between RNA methylation and DNA epigenetic mechanisms in plants [12]. However, there is a lack of knowledge regarding the crosstalk between 5mC and m 6 A regulators in HCC. Thus, the present work aimed to determine the crosstalk of 5mC and m 6 A in HCC, and develop an epigenetic and epitranscriptomic module eigengene (EME) reflecting 5mC and m 6 A modification levels that enabled to define clinical outcomes, therapeutic response and pharmacogenomic landscape.
The differential expression of 5mC/m 6 A regulators in tumors versus paired normal tissues; heterozygous and homozygous copy number variation (CNV), Pearson correlation between gene expression and CNV; the differential methylation in tumors versus paired normal tissues, the associations of methylation and expression, and prognosis influenced by methylation; spearman correlation analysis of the expression of 5mC/m 6 A regulators with pathway activity (activation or inhibition); and small molecule/drug sensitivity (half-maximal inhibitory concentration (IC50)) were analyzed via the Gene Set Cancer Analysis (GSCALite) web server based on the multi-omics data from 11160 samples across 33 TCGA cancer types, 746 drug sensitivity data from the Genomics of Drug Sensitivity in Cancer (GDSC) and the Cancer Therapeutics Response Portal (CTRP), and normal tissue expression data of 11688 samples from the GTEx [20].
Development of an epigenetic/epitranscriptomic module eigengene (EME) scoring system 5mC/m 6 A regulators were uploaded to the STRING online database (https:// string-db. org/) [21], and a protein-protein interaction (PPI) network was exported. To select hub 5mC/m 6 A regulators, the 41 5mC/m 6 A regulators were used as a module utilizing weighted gene co-expression network analysis (WGCNA) [22]. The summary expression level of the module was defined as the module eigengene via moduleEigengenes function. Then, the module membership (i.e., module eigengenebased intramodular connectivity) was computed as the association between the expression level of a given 5mC/ m 6 A regulator and the module eigengene. Hub 5mC/ m 6 A regulators were selected as those that had a module membership > 0.7. The overall expression value of the hub 5mC/m 6 A regulators was computed as the EME score. In accordance with the optimal cutoff of EME score, the patients were divided into low or high EME score group.

Functional enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) of hub 5mC/m 6 A regulators were analyzed utilizing clusterProfiler package [23]. The hallmark gene sets were acquired from the Molecular Signatures Database (http:// softw are. broad-insti tute. org/ gsea/ msigdb) [24]. Gene set variation analysis (GSVA) was adopted to ascribe signaling pathway variation scores to the gene sets, thus assessing the biological significance [25]. Gene Set Enrichment Analysis (GSEA) was also used to determine whether a defined gene set showed statistically significant between two groups [26].

Establishment of a nomogram
Independent risk factors derived from multivariate Cox regression analysis, were selected to establish a nomogram for predicting the likelihood of OS via rms package (version 6.2-0). Calibration curves were drawn to assess calibrating capacity. Area under the time-dependent receiver operating characteristic (ROC) curves (timedependent area under the curve (AUC)) and concordance index (C-index) were computed to assess discriminative capacity. The clinical benefits and utility of the nomogram were observed at different threshold probabilities by implementing decision curve analysis.

Therapeutic response analysis
Immunotherapy response was predicted via Immune Cell Abundance Identifier (ImmuCellAI) algorithm [31]. exclusion (TIDE) that may accurately estimate immunotherapy response [32]. The response to sorafenib (IC50) was predicted using pRRophetic package based on gene expression levels [33]. RNA-seq and clinical information of 58 patients who received anti-CTLA4 or anti-PD1 immunotherapy were collected from the GSE91061 cohort [34]. The predictive efficacy of the EME was verified in the cohort.

Terminal-deoxynucleoitidyl transferase mediated nick end labeling (TUNEL) staining
TUNEL staining was carried out utilizing TUNEL apoptosis kit (KTA2011; Abbkine) following the manufacturer's specifications. The percentage of TUNEL-positive cells was also calculated.

Flow cytometry for cell cycle analysis
Cells were seeded in six-well plates (1 × 10 5 cells / well). After 24 h, the cells were harvested after trypsinization and washed with ice-cold PBS. The cells were re-suspended in 300 μL PBS plus 5% bovine serum albumin (BSA; V900933; Sigma, USA) and fixed with 700 μL ethanol at 4 °C for 24 h after removing the supernatant following centrifugation. Next, the cells were washed with PBS and centrifuged to discard the ethanol, and re-suspended in 100 μL PBS plus 1 μL ribonuclease A (10 mg/ μL). After incubation for half an hour, the cells were stained with 50 μL propidium iodide at 37 °C for half an hour. Cell cycle distribution was evaluated with BD Biosciences FACSCalibur flow cytometry system and analyzed with FlowJo software (BD Biosciences, USA).

Wound healing
5 × 10 4 cells were planted onto 6-well plates. When the confluence was up to 90%, the cell monolayer was scratched with a 200 μL pipette tip. After removing detached cells, the remaining cells were cultivated in medium without FBS. Images were investigated at 0, and Development of an EME scoring system for reflecting 5mC/m 6 A modification status and predicting clinical outcomes in HCC. A Overview of the hub 5mC/m 6 A regulators and the EME scoring system. B The distribution of EME and clinical parameters across HCC patients. C The expression of the hub 5mC/m 6 A regulators along the EME. D PCA for the dissimilarity between the high and low EME subsets according to the expression matrix of the hub 5mC/m 6 A regulators. E Kaplan-Meier OS analysis for subgroup patients stratified by the EME. 24 h and then imaged. Cellular migration was quantified utilizing Image J software.

Transwell assays
For migration assay, 5 × 10 4 cells were suspended in 100 μL serum-free DMEM and planted onto the upper chamber of Transwell (

Single-cell sequencing (scRNA-seq) data
The scRNA-seq count matrix of 21 HCC samples was downloaded from the GSE149614 cohort [13]. Quality control was implemented with Seurat package [35].
Single cells with > 20% mitochondrial UMI counts as low-quality cells were removed. Batch effects were eliminated utilizing Integrate Data function. The top 15 principal components (PCs), and the top 1,500 highly variable genes, were selected. Influence of the percentage of mitochondrial UMI counts were removed with ScaleData function. Afterwards, cell populations were clustered with FindClusters function, and visualized with t-distributed stochastic neighbor embedding (t-SNE). The markers of each cell cluster were determined with FindAllMarkers function. The main cell types were determined on the basis of markers acquired from the CellMarker database (http:// biocc. hrbmu. edu. cn/ CellM arker/ or http:// bio-bigda ta. hrbmu. edu. cn/ CellM arker/) [36].

Cancer cell line data
We gathered drug sensitivity data of human cancer cell lines (CCLs) from the CTRP (https:// porta ls. broad insti tute. org/ ctrp) and PRISM (https:// depmap. org/ portal/ prism/) datasets that were acquired from the Cancer Cell Line Encyclopedia (CCLE) project (https:// porta ls. broad insti tute. org/ ccle/). Transcriptome data of the CCLE were employed for CTRP and PRISM analyses. The CERES score that can measure the dependency of genes in specific CCLs were downloaded from the dependency map (DepMap) portal (https:// depmap. org/ portal/), which was utilized for measuring the gene dependency in specific CCLs. Potential working mechanisms of the identified drugs were analyzed based on the Connectivity Map (CMap) database (https:// clue. io) [37].

Statistical analysis
Data processing, analysis, visualization, etc. were conducted with R packages (version 3.6.3; https:// www. bioco nduct or. org/) or Graphpad Prism software (version 8.0.1; https:// www. graph pad. com/). Differences between two groups were evaluated with student's t or Mann-Whitney U test. One-way analysis of variance (ANOVA) or Kruskal-Wallis test was conducted for multiple comparisons. Fisher's exact test was adopted for analyzing categorical data. Correlation between variables was determined through Pearson or Spearman test. Survival analysis was implemented by Kaplan-Meier approach and log-rank test through survival and survminer packages. Uni-and multivariable Cox regression models were built to calculate hazard ratios (HRs) and identify independent prognostic parameters. ROC curves were plotted using survivalROC package. Principal component analysis (PCA) was utilized to (See figure on next page.) Fig. 4 The excellent performance of the EME scoring system in predicting HCC progression. A Kaplan-Meier DSS analysis for subgroup patients stratified by the EME. B ROC curves of the EME for risk prediction of 1-, 3-and 5-year DSS. C Kaplan-Meier PFS analysis for subgroup patients stratified by the EME. D ROC curves of the EME for risk prediction of 1-, 3-and 5-year PFS. E-H Validation of OS outcomes of high and low EME subgroup patients and the prediction performance of the EME by ROC curves in the (E, F) LIRI-JP and (G, H) GSE15654 cohorts. I Distribution of the EME among normal liver, early-stage HCC, and advanced-stage HCC in the GSE6764 dataset. J-L ROC curves of the EME for distinguishing normal liver, early-stage HCC, and advanced-stage HCC in the GSE6764 dataset The EME scoring system accurately predicts the responses to mainstream therapies (TACE and sorafenib). A Comparison of the EME between subgroup patients stratified by clinicopathological parameters. B Difference in the EME between TACE responders and non-responders in the GSE104580 cohort. C ROC curves of the EME for evaluating the performance in predicting TACE response. D Difference in the EME between sorafenib responders and non-responders in the GSE109211 cohort. E ROCs of the EME for assessing the prediction value of the EME in sorafenib response. F Difference in the IC50 value of sorafenib between the high-and low-EME patients. G Heatmap of the expression of several sorafenib targets in the high-and low-EME patients. H, I The genomic alterations of sorafenib response-related genes in the high-and low-EME patients.
Percentage of mutations is presented on the right panel. Mutation burden is exhibited as a bar plot on the top. Mutation type is displayed in the bottom. *p-value < 0.05; ****p-value < 0.0001 assess the classification accuracy. P-value < 0.05 indicated statistical significance.

Multi-omics analysis of the crosstalk between 5mC and m 6 A regulators across pan-cancer
This study collected 21 and 20 genes that function as regulators of 5mC-epigenetic and m 6 A-epitranscriptomic mechanisms. The genome-wide omics data of the two regulator classes were analyzed in 33 cancer types. Firstly, we summarized the levels of 5mC and m 6 A regulators across pan-cancer. Most regulators exhibited the cooccurring genetic expression levels in each cancer type ( Supplementary Fig. 1A). Additionally, the expression levels of the two regulator classes were comparable across 33 cancer types. Through the GSCALite web server, we investigated the genetic alterations of 5mC and m 6 A regulators across pan-cancer. 5mC and m 6 A regulators presented comparable genetic alterations across 33 cancer types, and the remarkable co-occurrence of genetic alterations was found between the two regulator classes ( Supplementary Fig. 1B-E). We also investigated the DNA methylation features of 5mC and m 6 A regulators. Across most cancer types, there was a negative correlation between the mRNA expression and methylation for the same regulator (Supplementary Fig. 2A, B). Interestingly, the two regulator classes exhibited the comparable methylation levels across different cancer types, and their differential methylation status was correlated to survival outcomes ( Supplementary Fig. 2C-F). Genomic aberrations affect therapeutic responses and can be potentially applied for drug screening. We observed the significant correlations between drug sensitivity and mRNA expression of most 5mC and m 6 A regulators ( Supplementary  Fig. 3A, B).

Landscape of the crosstalk between 5mC and m 6 A regulators in HCC
We further evaluated the differential expression profiling of 5mC and m 6 A regulators in HCC using available tumor and normal tissue expression data. Almost all 5mC and m 6 A regulators were remarkably up-regulated in HCC than normal tissues ( Fig. 2A), which reflected the critical significance of epigenetic and epitranscriptomic regulation mechanisms in HCC. As illustrated in Fig. 2B and C, genetic mutations of 5mC and m 6 A regulators were not frequent in HCC. As for CNV characteristics, a majority of 5mC and m 6 A regulators displayed prevalent copy number amplifications or deletions (Fig. 2D, E). The above analysis revealed the transcriptomic and genomic characteristics of 5mC and m 6 A regulators in HCC. Further investigation showed that the 5mC and m 6 A regulators interacted with one another frequently (Fig. 2F). Moreover, at the transcriptomic levels, positive interactions between 5mC and m 6 A regulators were observed (Fig. 2G). For example, DNMT1 and METTL3 were both highly expressed in two HCC cell lines (HepG2 and Huh-7) than normal hepatocytes L-02 (Fig. 2H, I). Silencing DNMT1 or METTL3 significantly attenuated the expression of the other in HepG2 and Huh-7 cell lines (Fig. 2J-O). Similarly, it was proven that the expression of TET3 and YTHDC1 influenced each other in HCC cells (Supplementary Fig. 4A-F).

Development of an EME scoring system for reflecting 5mC/ m 6 A modification status in HCC
Pan-cancer survival analysis showed that 5mC and m 6 A regulators played similar roles in OS outcomes of pan-cancer ( Supplementary Fig. 5A). Most regulators served as risk factors of HCC patients' OS (Supplementary Fig. 5B). To determine the hub regulators involved in 5mC and m 6 A modifications in HCC, we employed WGCNA to select hub 5mC/m 6 A regulators that exhibited the remarkably high correlations, which can be explained by their crosstalk. Functional enrichment analysis proved the key roles of hub 5mC/m 6 A regulators in epigenetic and epitranscriptomic processes ( Supplementary Fig. 6A-D). An EME scoring system was developed to reflect both the 5mC-epigenetic and m 6 A-epitranscriptomic modification levels through calculating the overall expression levels of hub 5mC/m 6 A regulators (Fig. 3A). According to the optimal cutoff value (0.004) of EMEs, TCGA-LIHC cohort was separated into the high-and low-EME subsets (Fig. 3B, C). The t-SNE highlighted the distinct transcriptional program for two subsets defined by hub 5mC/m 6 A regulators (Fig. 3D).

The EME scoring system for risk stratification and OS outcomes of individual HCC patients
The clinical outcomes of high-and low-EME HCC patients remarkably varied. Patient with high EME had shorter OS relative to those with low EME (Fig. 3E). To quantify the capacity of this scoring system to predict OS, ROC curves were generated. The 1-, 3-, and 5-year AUC values were 0.71, 0.63 and 0.64, respectively, implying that the scoring system performed well and had the excellent predictive ability (Fig. 3F). To examine the independence of the EME, uni-and multivariate cox regression analyses were implemented. After adjusting by common clinical parameters, the EME still exhibited a robust assessment capacity in OS (Fig. 3G). To further confirm the robustness and superiority of the EME, we implemented subgroup analysis according to age, gender, grade, stage, T stage, and HBV infection. In all subgroups, patients with high EME presented shorter OS than those with low EME (Supplementary Fig. 7A-L). Overall, the EME could greatly optimize risk stratification as well as accurately predict HCC prognosis.
To quantify the risk evaluation for individual HCC patients, a personalized scoring nomogram that combined the EME with clinicopathological parameters (age and stage) was generated for predicting 1-, 3-and 5-year OS probability (Fig. 3H). Based upon the calibration curves, the nomogram exhibited the excellent performance in comparison to an ideal model (Fig. 3I). ROC curves also showed that the nomogram had the strong clinical significance for predicting OS (Fig. 3J). Moreover, the C-index was 0.668. In summary, the nomogram for OS had a considerable discriminative and calibrating capacity. The clinical benefits of the nomogram were also evaluated by decision curve analysis. In Fig. 3K, the nomogram enabled to better predict 1-, 3-, and 5-year OS, because it added more net benefits for almost all threshold probabilities.

The excellent performance of the EME scoring system in predicting HCC progression
Based on the robust association of the EME with OS, we hypothesized that it would have the same relationships with disease-specific survival (DSS) and PFS. We found that the patients with low EME had a remarkable survival advantage, and the EME enabled to predict the DSS and PFS time (Fig. 4A-D). Aiming to verify the prediction value of the EME, we enrolled the LIRI-JP and GSE15654 independent cohorts. High-EME patients' OS was remarkably better; additionally, the EME exhibited the excellent predictive power for OS (Fig. 4E-H). The GSE6764 dataset was employed to evaluate whether the EME was associated with the pathological progression of HCC. As the disease progressed, the EME increased gradually (Fig. 4I). ROC curves demonstrated that the EME could accurately differentiate distinct pathological stages of HCC (Fig. 4J-L). Higher EME was found to be linked to advanced pathological grade and stage (Fig. 5A). Altogether, the EME can reflect the progression of HCC.

The EME scoring system accurately predicts the responses to mainstream therapies (TACE and sorafenib)
In addition to surgical resection, the mainstream therapies against HCC mainly comprise TACE and sorafenib. We evaluated the capacity of the EME scoring system to predict the responses to TACE and sorafenib. In the GSE104580 cohort, TACE non-responders exhibited the higher EME than TACE responders (Fig. 5B). ROC curves were generated to reflect the predictive accuracy. The EME exhibited the high AUC value of 0.722 (Fig. 5C), implying that the scoring system possessed the excellent performance in predicting the response to TACE. Next, we calculated the difference in the EME between sorafenib responders and non-responders in the GSE109211 dataset. As shown in Fig. 5D, we found that the EME score of sorafenib non-responders was higher than that of sorafenib responders. The AUC value was 0.894 (Fig. 5E), demonstrating that the EME had the well accuracy in predicting the response to sorafenib. Further analysis using pRRophetic package showed that the low-EME patients were more sensitive to sorafenib (Fig. 5F). Additionally, we found that several sorafenib targets (BRAF, RAF1, FLT1, KIT, FLT3, FLT4) had the higher expression in the high-EME patients (Fig. 5G). The higher mutation burden of sorafenib response-related genes was found in the low-EME patients (Fig. 5H, I).

5mC and m 6 A regulators cooperatively weaken apoptosis and facilitate proliferation for HCC cells
Based on the relevance of the EME scoring system to HCC progression, we next explored the mechanisms by which the interplay between 5mC and m 6 A regulators regulated HCC progression. We found that most 5mC and m 6 A regulators were significantly correlated to tumorigenic pathways (Fig. 6A). Furthermore, the EME score was positively linked with tumorigenic pathways (such as epithelial-to-mesenchymal transition (EMT), DNA damage, and apoptosis; Fig. 6B). Herein, we focused on the synergistic effect of 5mC regulator DNMT1 and m 6 A regulator METTL3 on apoptosis and proliferation of HCC cells. TUNEL staining showed that apoptotic level of HCC cells was enhanced by si-DNMT1 or si-METTL3 in HCC cells (Fig. 6C-E). Simultaneous inhibition of DNMT1 and METTL3 synergistically enhanced apoptosis in HCC cells. We also conducted colony formation assay and EdU staining to evaluate the proliferative capacity of HCC cells. It was found that si-DNMT1 or si-METTL3 mitigated the proliferative capacity of HepG2 and Huh-7 cells, and synergistic effect was observed when DNMT1 and METTL3 were simultaneously inhibited ( Fig. 6F-K). Above evidence implied that 5mC and m 6 A regulators cooperatively weakened apoptosis and facilitated proliferation for HCC cells.

5mC and m 6 A regulators cooperatively enhance DNA damage repair and cell cycle progression in HCC
DNA damage response is an important signaling process regulating DNA repair, cell cycle arrest, and maintaining genome homeostasis in response to endogenous and exogenous genotoxic stress. Evidence suggests that HCC cells display a strong DNA repair ability to overcome DNA damage caused by treatment [38]. In Fig. 7A, B, we found that DNA damage repair, cell cycle-related pathways (G2M checkpoint, mitotic spindle, E2F, P53, MYC, PI3K Akt mTOR signaling, TGF-β signaling, etc.) exhibited the strong activation in the high EME patients, implying the cooperative effect of 5mC and m 6 A regulators on DNA damage repair and cell cycle progression. For this reason, we further analyzed the associations between the EME and different DNA damage repair signatures. Base excision repair, nonhomologous end joining, homologous recombination, mismatch repair, and Fanconi anemia were positively correlated to the EME (Fig. 7C, D), and these pathways were significantly enriched in the high EME group (Fig. 7E), implying that the crosstalk between 5mC and m 6 A regulators might be functionally essential for DNA damage response. We examined the effects of knockdown of DNMT1 and METTL3 on cell cycle through flow cytometry. Compared with knockdown of DNMT1 or METTL3 alone, we found that simultaneous inhibition of DNMT1 and METTL3 synergistically induced G2/M cell cycle arrest in HCC cells (Fig. 7F-H).

5mC and m 6 A regulators synergistically accelerate migration, invasion and EMT of HCC cells
Further analysis showed that the EME was positively linked to EMT process (Fig. 8A-C). Wound healing and transwell assays showed that simultaneous inhibition of DNMT1 and METTL3 synergistically attenuated migration and invasion in HCC cells (Fig. 8D-L). Additionally, simultaneous inhibition of DNMT1 and METTL3 synergistically lowered the expression of β-catenin and enhanced the expression of E-cadherin both in HepG2 and Huh-7 cells (Fig. 8M-R), implying the inactivation of EMT process. In summary, 5mC and m 6 A regulators might synergistically accelerate HCC metastasis.

Correlation between the EME scoring system and potential extrinsic and intrinsic immune escape mechanisms in HCC
Immune surveillance is a key mechanism to prevent tumor development and progression. Extrinsic factors (immune cells) and intrinsic factors (including gene alterations and pathway activity) play important roles in immune escape mechanisms. Hub 5mC/m 6 A regulators were remarkably linked to most tumor-infiltrating immune cells (Fig. 9A). ScRNA-seq can maximize unbiased information to explore transcriptional diversity at the single-cell level. A total of 21 HCC scRNA samples were involved in our study. After quality control, normalization, and dimensionality reduction analysis (Supplementary Fig. 8A-G), the cells were classified as 9 main cell lineages (Fig. 9B), comprising NK cells, monocytes, iPS cells, T cells, hepatocytes, endothelial cells, smooth muscle cells, macrophages, and B cells. Figure 9C visualized the top five markers of the main cell lineages. We found that 5mC/m 6 A regulators were almost expressed in all immune and non-immune cell lineages (Fig. 9D, E). Through implementing different algorithms, we estimated the infiltration levels of immune cells across HCC. As a result, the EME was remarkably linked to the infiltration levels of most immune cell types (Fig. 9F-H), implying the widespread 5mC and m 6 A modifications. Both in TCGA-LIHC and LIRI-JP datasets (Fig. 9I, J), the EME presented positive correlations to macrophage markers (such as FIZ1, TGFB1, IL15RA, and IL12A). In addition, a majority of MHC molecules and immune checkpoints (Fig. 9K) as well as immunomodulators (chemokines, receptors, MHC, and immuno-stimulators; Fig. 9L) were positively correlated to the EME. CD8 + T cell exhaustion is a primary barrier to current immunotherapy [39]. The EME was positively correlated to CD8 + T cell exhaustion markers (CD276, TIGIT, LAG3, CD27, CXCL9, IDO1, etc.; Fig. 9M). Additionally, the remarkable difference in the EME was found among known immune subtypes, with the highest EME in C1 (Fig. 9N). Persistent expression of PD-L1 in tumor cells promotes tumor cells escape from immune surveillance and host T cell exhaustion [40]. Both si-DNMT1 and si-METTL3 mitigated the expression of PD-L1 in HepG2 and Huh-7 cells (Fig. 9O-Q). Simultaneous inhibition of DNMT1 and METTL3 synergistically lowered the expression of PD-L1. In summary, EME scoring system was remarkably linked to extrinsic and intrinsic immune escape mechanisms in HCC.

The EME scoring system accurately predicts immunotherapeutic response and hyper-progression in HCC
Based on the aforementioned strong relationship between the EME and antitumor immunity, we further evaluated whether the EME could be used to predict the response to immunotherapy. The high-EME subset had the higher enrichment scores of most immunotherapypredicted gene signatures than the low-EME subset (Fig. 10A), which was confirmed in the LIRI-JP cohort (Fig. 10B). Moreover, we found that the EME exhibited the positive correlations to most immunotherapy-predicted gene signatures (Fig. 10C). The cancer immune cycle includes a series of events required for immunemediated tumor growth control. Disruption of one or more steps enables tumor cells to escape immune surveillance. The EME was significantly linked to most steps within cancer immune cycle (Fig. 10C). Additionally, the high-EME subset had the higher enrichment scores of immune inhibited oncogenic pathways and EGFR ligands (Fig. 10D). We also computed the TIDE score across HCC, a reliable immunotherapy response predictor. The higher TIDE was observed in the high-EME subset (Fig. 10E), and the TIDE was positively correlated to the EME (Fig. 10F), implying that the low-EME patients were more likely responsive to immunotherapy. Immu-CellAI was also employed to predict cancer immunotherapy response. Immunotherapy non-responders had the higher EME than responders (Fig. 10G). The AUC value was 0.634, implying the well predictive accuracy of the EME in immunotherapy response (Fig. 10H). Immunotherapy response was also predicted by TIDE algorithm. The higher EME was observed in immunotherapy non-responders than responders both in TCGA-LIHC and ICGC cohorts (Fig. 10I, J). The efficacy of the EME in predicting immunotherapy response was also validated in an anti-CTLA4 and anti-PD1 immunotherapy cohort (GSE91061). We computed the EME of patients in this cohort, and assessed the difference in the EME among patients who differently responded to immunotherapy. The low-EME subset presented a higher proportion of complete response (CR) / partial response (PR) compared with the high EME subset (24% vs. 15%; Fig. 10K), and patients who completely responded to anti-CTLA4 or anti-PD1 therapy exhibited the lowest EME (Fig. 10L). Additionally, the low-EME patients had a remarkable survival advantage (Fig. 10M). These results implied that the EME scoring system enabled to reflect the patients' sensitivity to immunotherapy. We also focused on the immunotherapy-associated hyperprogression. The high-EME subset displayed the higher transcript levels and copy number amplification frequencies of genes positively linked to hyper-progression, such as CCND1, EGFR, and FGF4 (Fig. 10N, O). In summary, the high-EME HCC patients might not benefit from immunotherapy and instead present a higher possibility of hyper-progression.

Correlation of genomic alterations with the EME in HCC
The relationship between the EME and genomic alterations was further analyzed. Firstly, we observed that the high-EME subset had the higher levels of immunogenicity indicators (aneuploidy score, CTA score, HRD score, and intratumor heterogeneity; Fig. 11A-D). Next, we computed the GISTIC score in the high-and low-EME patients. Compared with the high-EME subset, higher GISTIC score and CNV frequency were found in the low-EME subset (Fig. 11E, F). Then, FGA, FGG, and FGL were calculated, respectively. The low-EME subset displayed higher FGA, FGG and FGL levels than the high-EME subset (Fig. 11G, H). This implied that the high 5mC/m 6 A modification might contribute to a reduced copy number gain/loss frequency in HCC. Next, we analyzed the difference in gene mutations between the highand low-EME subsets. However, we found that there was no remarkable difference between the two subsets (Fig. 11I).

Identification of potential therapeutic compounds and druggable targets for HCC patients with high EME
The drug response datasets CTRP and PRISM were employed to select potential compounds for the high-EME subset. The compounds with lower AUC in this subset were screened via evaluating the differences in drug response between the high-and low-EME subsets. As a result, five CTRP-derived compounds (paclitaxel, austocystin D, cucurbitacin I, GSK461364, and SB-743921; Fig. 11J) and six PRISM-derived compounds (norfloxacin, LY2606368, mitoxantrone, volasertib, navitoclax, and sirolimus; Fig. 11K) were determined via Spearman correlation analysis between the EME and AUC value (p-value < 0.05 and Spearman's r < − 0.3). In Fig. 11L, we showed the working mechanisms of above compounds: Tubulin inhibitor for paclitaxel; inducer of DNA damage for austocystin D; JAK inhibitor for cucurbitacin I; PLK inhibitor for GSK461364; KIF11 inhibitor for SB-743921; bacterial enzyme DNA gyrase inhibitor for norfloxacin; CHK1 inhibitor for LY2606368; topoisomerase inhibitor for mitoxantrone; PLK1 inhibitor for volasertib; BCL inhibitor for navitoclax; and mTOR inhibitor for sirolimus. Proteins that were highly positively linked to the EME might exhibit potential therapeutic implications for the high-EME patients. Nevertheless, most human proteins are still undruggable because they are short of distinct active sites to which compounds are capable of binding or residing within cells that are inaccessible to biological agents.
Thus, to determine potential druggable targets for the high-EME patients with undesirable survival outcomes, this study obtained the target information of 6,125 compounds. Firstly, we computed Spearman's correlation of the protein levels of druggable targets with the EME. As a result, 671 druggable targets were determined (p-value < 0.05 and Spearman's r > 0.4; Fig. 11M). Then, by implementing correlation analyses on the CERES score of druggable targets and the EME, we determined 497 druggable targets (p-value < 0.05 and Spearman's r < -0.45; Fig. 11N). Based on the previously identified druggable targets of HCC [41], two druggable targets including TUBB1 and P2RY4, were finally determined by above analysis. In HCC cells, silencing DNMT1 or METTL3 significantly decreased the expression of TUBB1 and P2RY4 (Fig. 11O-S). Simultaneous suppression of DNMT1 and METTL3 synergistically reduced their expression, implying that suppressing the function of these two targets in the high-EME subset might contribute to beneficial therapeutic response.

Discussion
The present study unveiled the crosstalk between 5mCepigenetic mechanism and m 6 A-epitranscriptomic mechanism in HCC at the multiomic levels. The EME scoring system was developed to reflect 5mC/m 6 A modification status in HCC, which may greatly optimize risk stratification, and accurately predict HCC patients' responses to mainstream therapies (TACE and sorafenib) and immunotherapy as well as hyper-progression. Moreover, the crosstalk between 5mC and m 6 A modifications defined pharmacogenomic landscape including potential therapeutic compounds and druggable targets (TUBB1 and P2RY4). Altogether, our findings laid a solid foundation for epigenetic regulation of HCC as well as paved a new avenue for relevant therapeutic targets.
Interference with tumor growth via 5mC or m 6 A modification levels is a promising therapeutic strategy for HCC. Phenotypic intra-tumor heterogeneity explains the poor efficacy of single-target systemic therapies in HCC [42]. Suppression of 5mC regulator (DNMT1) and m 6 A (See figure on next page.) Fig. 10 The EME scoring system accurately predicts immunotherapeutic response and hyper-progression in HCC. A, B Comparison of the enrichment scores of immunotherapy-predicted signatures in the high-and low-EME patients in TCGA-LIHC and LIRI-JP cohorts. C Correlations between the EME and immunotherapy-predicted pathways (left panel) and cancer immune cycle (right panel). D Comparison of the activity of immune inhibited oncogenic signaling and EGFR ligands in the high-and low-EME subsets. E Difference in the TIDE score between the high-and low-EME subsets. F Correlation between the EME and TIDE score. G Difference in the EME between ImmuCellAI-predicted immunotherapy responders and non-responders. H ROC curves for evaluating the predictive accuracy of the EME in immunotherapy response. I, J Comparison of the EME between TIDE-predicted immunotherapy responders and non-responders in TCGA-LIHC and ICGC cohorts. K The proportion of patients in the GSE91061 cohort with different clinical responses (CR, complete response; PR, partial response; SD, stable disease; and PD, progressed disease) in the high-and low-EME subsets. L Comparison of the EME among patients with different clinical responses in the GSE91061 cohort. M Kaplan-Meier curves of OS between the low-and high-EME patients in the GSE91061 cohort. N, O Comparison of the mRNA expression and copy number amplification/deletion frequencies of hyper-progression-associated genes in the low-and high-EME patients. Ns: no significance; ***p-value < 0.001; ****p-value < 0. regulator (METTL3) cooperatively attenuated HCC progression. HCC tumors exhibit intra-and intertumoral heterogeneity at the molecular, histological, and clinical levels. As a multifocal neoplasm, there are differences between distinct lesions in the same patient, known as "interlesion" heterogeneity. The complex heterogeneity hinders the development of HCC treatment, notably those who respond differently to drugs and immunotherapies. Hence, stratification of HCC tumors into clinically and molecularly homogeneous subgroups may improve physicians' therapeutic options. In the present study, the simple and easily applicable EME scoring system was developed to predict clinical outcomes as well as treatment responses (TACE and sorafenib) in patients with HCC, which was verified in the independent and external cohorts.
To better understand the cell cycle, the central molecules that participate in this process have been discovered, along with the introduction of the checkpoint concept [43]. The regulatory network comprising the central molecules is capable of accurately regulating the process. The shared biological feature of HCC is uncontrolled growth, manifested by disordered cell cycle resulting in uncontrollable cell proliferation and attenuated apoptosis [38]. DNA damage repair is a crucial event in mediating cell cycle arrest. The work demonstrated that the EME exhibited positive correlations to diverse DNA damage repair mechanisms, and suppression of DNMT1 and METTL3 synergistically mitigated cell cycle progression of HCC cells.
Complex genomics and tumor microenvironment create a molecular conundrum for the diagnosis and therapy of HCC and result in therapeutic failure and eventually fatal outcomes [44]. A previous study classified liver cancer into four immune subtypes: tumorassociated macrophages, CTNNB1, cytolytic activity (CYT), and regulatory T cells (Tregs) [45]. Among them, CYT and Tregs subtypes are inflamed tumors, and TAMs and CTNNB1 subtypes are non-inflamed tumors. Emerging scRNA-seq can obtain quantitative transcriptomic data at the single-cell level, while eliminating the bias that occur in RNA-seq caused by diverse cell subpopulations [46]. The 5mC and m 6 A regulators were extensively distributed in most immune and nonimmune cell populations across HCC. Although HCC tumors harbor significant infiltration of immune cells, they cannot kill tumor cells by inference [47]. Epigenetic and epitranscriptomic modifications can impact the functions and phenotypes of immune cell types for cell killing and functional tuning. The transcriptional and epigenetic landscape of exhausted CD8 + T cells defines exhaustion as a distinct branch of CD8 + T cell differentiation. Selective epigenetic reprogramming alters the T-cell landscape in HCC and may enhance the therapeutic efficacy. T cells cannot infiltrate tumors due to vasculature triggered by tumor cells, chemokines as well as known immunosuppressive cells (MDSCs, immature DCs, macrophages, Tregs, etc.). Epigenetic modulation involves all three aspects. PD-L1 binds to its receptor PD-1, thus leading to immune escape by counteraction of activating signaling on T cells [48]. Moreover, CD80 can bind to PD-L1 and transmit negative signaling. PD-L1 is upregulated in response to some inflammatory signals (IFN-γ, etc.) generated by active T cells during antitumor immune responses [49]. Simultaneous inhibition of DNMT1 and METTL3 synergistically lowered the expression of PD-L1 in HCC cells, thereby suppressing immune escape. Immunosuppression and immune evasion exert key roles in tumorigenesis and tumor progression, in which tumor cells can innately or adaptively express immunosuppressive molecules to escape host immune attack [50]. An in-depth understanding of immune escape mechanisms in tumor cells is critical for overcoming resistance as well as enabling innovative progress in immunotherapy. The advent of ICIs has notably altered the landscape of HCC therapy. Despite these improvements, one of the most relevant unaddressed medical requirements in the field is to identify a biomarker of therapeutic response that may assist determine patients who can respond to ICIs. The EME may predict efficacy of ICI therapy for HCC. From a methodological point of view, the EME scoring system may assist to identify potential responders to immunotherapy and reduce side effects for patients who are not likely to benefit from it. Fig. 11 Associations between the EME and genomic alterations, potential therapeutic compounds and druggable targets in HCC. A-D Comparison of aneuploidy score, CTA score, HRD score, and intratumor heterogeneity in the low-and high-EME patients. E, F CNVs in the low-and high-EME patients. G, H Differences in FGA, FGG, and FGL between the high-and low-EME patients. I The top twenty mutated genes across HCC patients stratified by the EME. J, K Spearman correlation analysis on the EME with AUC of CTRP-and PRISM-derived compounds (left), and comparison of AUC between the low-and high-EME subsets (right). L Potential working mechanisms of the identified compounds. M Spearman's correlation on the protein expression of compound targets and the EME. Blue dot represents a significant positive correlation (p-value < 0.05 and Spearman's r > 0.4). N Spearman's correlation on the CERES score of druggable targets and the EME. Red dot denotes a significant negative correlation (p-value < 0.05 and Spearman's r < -0.45). O-S Western blot for the expression of TUBB1 and P2RY4 in the presence or absence of si-DNMT1 or si-METTL3 in HepG2 and Huh-7 cells. Ns: no significance; *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001